SedimentRouting Subroutine

public subroutine SedimentRouting(dt, flow, v, dp)

compute sediment routing in channelized drainage network using Muskingum-Cunge method. Interrill erosion is supposed to act as suspended sediment so deposition is not contemplated. Bedload sediment is divided in different grainsize classes and routing is performed for each class separately. Deposition is contemplated comparing sediment discharge to transport capacity. Flow Routing method code remind: 0 --> hillslope, water flow rate is not defined 1,2,3,11,30 --> channel routing, water flow is computed according to different schemes 5 --> Instantaneous mass transferring inside lakes (infinitive celerity) 1000:2000 --> reservoir

References:

Sun, H., P. S. Cornish, and T. M. Daniell, Contour-based digital elevation modeling of watershed erosion and sedimentation: Erosion and sedimentation estimation tool (EROSET), Water Resour. Res., 38(11), 1233, doi:10.1029/2001WR000960, 2002.

Singh, V.P., Quiroga, C.A., A dam-breach erosion model: I. formulation, Water resources management, 1, 177-197, 1987.

Arguments

Type IntentOptional Attributes Name
integer(kind=short) :: dt

routing time step
type(grid_real), intent(in) :: flow
type(grid_real), intent(in) :: v

water flow velocity in channel (m/s)
type(grid_real), intent(in) :: dp

water depth (m)

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: cappa
real(kind=float), public :: ddx
integer(kind=short), public :: iin
integer(kind=short), public :: is
integer(kind=short), public :: jin
integer(kind=short), public :: js
integer(kind=short), public :: k
real(kind=float), public :: transportCapacity
real(kind=float), public :: vel
real(kind=float), public :: wd
real(kind=float), public :: width

Source Code

SUBROUTINE SedimentRouting &
!
(dt, flow, v, dp)


!Arguments with intent in:
INTEGER (KIND = short) :: dt !!routing time step
TYPE (grid_real), INTENT(IN) :: flow
TYPE (grid_real), INTENT(IN) :: v !!water flow velocity in channel (m/s)
TYPE (grid_real), INTENT(IN) :: dp !!water depth (m)

!local declarations:
INTEGER (KIND = short) :: k, iin, jin, is, js
REAL (KIND = float) :: ddx, cappa
REAL (KIND = float) :: vel
REAL (KIND = float) :: transportCapacity
REAL (KIND = float) :: wd
REAL (KIND = float) :: width

!------------end of declaration------------------------------------------------

QinSS = 0.

DO K=1,sedReach%nreach
  iin=sedReach%branch(k)%i0
  jin=sedReach%branch(k)%j0
  DO WHILE ( .NOT.((jin == sedReach%branch(k)%j1) .AND. &
	               (iin == sedReach%branch(k)%i1)) )
	               
     CALL DownstreamCell (iin, jin, sedFlowDirection % mat(iin,jin), &
                          is, js, ddx, sedFlowDirection)
     
     SELECT CASE (sedReach % branch (k) % routingMethod)
     
       CASE (0)
          !compute cumulated sediment rate on hillslope along drainage network 
          !first cell of order 1 reach
          IF (sedReach % branch (k) % order == 1 .AND. &
              sedReach % branch (k) % I0 == iin   .AND. &
              sedReach % branch (k) % J0 == jin) THEN
             QinSS % mat(iin,jin) = interrillErosion % mat(iin,jin)
             QinSS % mat(is,js) =  QinSS % mat(iin,jin) + interrillErosion % mat(is,js)
          ELSE
             QinSS % mat(is,js) = QinSS % mat(is,js) + QinSS % mat(iin,jin) +  interrillErosion % mat(is,js)
          END IF
		  !QinSS % mat(is,js) = QinSS % mat(is,js) + QinSS % mat(iin,jin)
		  
		  !update storage sediment variation. On hillslope variation is always negative
		  deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) - &
		                             interrillErosion % mat(iin,jin) * dt / 1000000.
		                             

       CASE (1,2,3,11,30)
          ! sediment routing in channel reach  
          
          !suspended solid
          !^^^^^^^^^^^^^^^
!          QoutSS % mat(iin,jin) = QinSS % mat(iin,jin)  * gridC1 % mat(iin,jin) + &
!                                  PinSS % mat(iin,jin)  * gridC2 % mat(iin,jin) + &
!                                  PoutSS % mat(iin,jin) * gridC3 % mat(iin,jin)
      
           !DA SISTEMARE - GR 12/10/2016
           !vel = flow % mat (iin,jin) / (dp % mat (iin,jin) * dp % mat (iin,jin) * sedReach % branch (k) % B0) 
          
           IF (vel > 0.) THEN
 
              cappa = ddx / vel
              !cappa = 3
              
              CALL MuskingumCungeTodini ( dt = dt, dx = ddx, &
                                 Qin =  QinSS%mat(iin,jin), &
                                 Pin =  PinSS%mat(iin,jin), &
                                 Pout =  PoutSS%mat(iin,jin), &
                                 Qlat = 0., Plat = 0., &
                                 B0 = 0., alpha = 0., slope = 0., n = 0., &
                                 Qout = QoutSS%mat(iin,jin), &
                                 depth = wd, width =  width, k = cappa, x = 0. )
       
  
           ELSE
              QoutSS % mat(iin,jin) = 0.
           END IF

          !filter negative discharge                                            
          IF ( QoutSS % mat(iin,jin) < 0.) THEN
             QoutSS % mat(iin,jin) = 0.
          END IF
          
           !limit suspended load to transport capacity
           !assume particle size coming from hillslope of 0.02 mm
          transportCapacity = Yang1979 (flow % mat(iin,jin),       &
                                        sedReach % branch(k) % slope, &
                                        vel, 0.02, &
                                        dp % mat(iin,jin)          )
          IF ( QoutSS % mat(iin,jin) > transportCapacity ) THEN
            
!             deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) + &
!		                            (QoutSS % mat(iin,jin) - transportCapacity) * dt / 1000000.
             QoutSS % mat(iin,jin) = transportCapacity
             QinSS % mat(iin,jin) =  transportCapacity
             
          END IF
          
           !update storage sediment variation
          deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) + &
		                            (QinSS % mat(iin,jin) - QoutSS % mat(iin,jin)) * dt / 1000000.
          
          
          !update downstream input discharge
          QinSS % mat(is,js) =  QinSS % mat(is,js) + QoutSS % mat(iin,jin)
          !store computed discharge for next time step
          PinSS % mat(iin,jin) = QinSS % mat(iin,jin)
          PoutSS % mat(iin,jin) = QoutSS % mat(iin,jin)
          
   
          !bed load transport
          !^^^^^^^^^^^^^^^^^^
          !compute source term
          QinBL % mat(iin,jin) = QinBL % mat(iin,jin) + ChannelDetachmentRate (flow % mat(iin,jin), &
                                                  sedReach%branch(k)%slope, &
                                                  rusleK % mat (iin,jin),    &
                                                  rusleC % mat (iin,jin)     )
                                                  
          !update storage sediment variation
!		  deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) - &
!		                             QinBL % mat(iin,jin) * dt / 1000000.
                                                                 
          !route downstream                                        
!          QoutBL % mat(iin,jin) = QinBL % mat(iin,jin)  * gridC1 % mat(iin,jin) + &
!                              PinBL % mat(iin,jin)  * gridC2 % mat(iin,jin) + &
!                              PoutBL % mat(iin,jin) * gridC3 % mat(iin,jin)
                              
            IF (vel > 0.) THEN
              cappa = ddx / vel
              
              CALL MuskingumCungeTodini ( dt = dt, dx = ddx, &
                                 Qin =  QinBL%mat(iin,jin), &
                                 Pin =  PinBL%mat(iin,jin), &
                                 Pout =  PoutBL%mat(iin,jin), &
                                 Qlat = 0., Plat = 0., &
                                 B0 = 0., alpha = 0., slope = 0., n = 0., &
                                 Qout = QoutBL%mat(iin,jin), &
                                 depth = wd, width = width, k = cappa, x = 0. )
  
           ELSE
              QoutBL % mat(iin,jin) = 0.
           END IF                   
                              
                              
                              
          !filter negative discharge                                            
          IF ( QoutBL % mat(iin,jin) < 0.) THEN
             QoutBL % mat(iin,jin) = 0.
          END IF 
                                
          !limit bed load to transport capacity
          transportCapacity = Yang1973 (flow % mat(iin,jin),       &
                                        sedReach%branch(k)%slope, &
                                        vel, sedReach%branch(k)%d50, &
                                        dp % mat(iin,jin)          )
          IF ( QoutBL % mat(iin,jin) > transportCapacity ) THEN
             !update storage sediment variation if deposition is positive
!             deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) + &
!                     (QoutBL % mat(iin,jin) - transportCapacity) * dt / 1000000.
             QoutBL % mat(iin,jin) = transportCapacity
             QinBL % mat(iin,jin) = transportCapacity
             
          END IF
          !update storage sediment variation 
           deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) + &
		                            (QinBL % mat(iin,jin) - QoutBL % mat(iin,jin)) * dt / 1000000.
          
                                                            
          !update downstream input discharge
          QinBL % mat(is,js) =  QinBL % mat(is,js) + QoutBL % mat(iin,jin)
          !store computed discharge for next time step
          PinBL % mat(iin,jin) = QinBL % mat(iin,jin)
          PoutBL % mat(iin,jin) = QoutBL % mat(iin,jin)
          
          !build QoutSed
          !^^^^^^^^^^^^^
          QoutSed  % mat(iin,jin) = QoutSS % mat(iin,jin)  + QoutBL % mat(iin,jin)
          !QoutSed  % mat(iin,jin) = QoutBL % mat(iin,jin)
          
          
     END SELECT
     
     !go downstream 
     jin = js
     iin = is
  END DO
  !last cell of last reach
  IF (K == sedReach%nreach) THEN
      
  END IF
END DO


END SUBROUTINE SedimentRouting